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Abstract: The conventional approximate formula for neutrino oscillation in matter 

which is obtained from the expansion in terms of the ratio of mass square differences 
a = Am 2 i/Am|^ ~ 0.03, first proposed by Cervera, et al and Freund, turns out to be an 
accurate formula for accelerator neutrino experiments. Originally it required the neutrino 
energy to be well above the solar resonance to validate the expansion but it is found to be 
still very accurate when the formula is extrapolated to the resonance, which is practically 
important for the T2K experiment. This paper shows that the accuracy is guaranteed by 
cancellations of branch cut singularities and also, for the first time, analytically computes 
the actual error of the formula. The actual error implies that the original requirement can 
be safely removed in current experiments. 
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1 Introduction 


In long-baseline(LBL) neutrino experiments, the matter effect [1-3] is usually not negligi¬ 
ble. For current LBL accelerator neutrino experiments such as T2K[4, 5], MINOS[6] and 
NOvA[7, 8] where the matter densities are almost constant, there is a useful approximate 
formula for the transition probability. Taking the same notations as PDG, the formula is[9] 
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( 1 . 1 ) 


where 

A = 2\/2GirAeP/Amii, a = Aml^/Aml^ ^ 0.03, (1.2) 

and A = Am^iL/(4E). is the electron number density in matter, about 1.4cm“^ACi in 
the Earth’s crust. 

The formula was originally derived in [10, 11] as a series expansion in a. But the 
problem is that due to the non-perturbative behavior near the solar resonance, the expansion 
is expected to be valid only when the neutrino energy is well above the solar resonance. 
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This was emphasized in ref. [11], because the approximation a/A<^\ was used when the 
formula was derived. We will reformulate the derivation of the formula in section 2 to show 
the problem more explicitly but here we take the solar mixing angle 9i2 as a good example 
to show the problem. The effective sin20i2 in matter, denoted as sin20]^, expanded in a 
to first order, is [11] 

sin20[^~a/A. (1.4) 

The solar resonance is at ^ = a cos 20i2 ~ 0.4a so near the solar resonance sin 20]^ is quite 
likely to be larger than 1. As will be shown in section 2, sin 20]^ > 1 does appear in the 
expansion when the energy is lower than 0.34GeV, which makes the calculation invalid. 
Originally, sin 20]^ in the calculation was expected not only less than 1 but also small, i.e. 
sin20[^ <C 1, otherwise the unitarity of the effective mixing matrix will be badly violated, 
thereby invalidating the calculation. 

Despite the claimed bound (1.3) in [11], in practice this formula works well below the 
bound (see figures 6, 7 presented later in this paper). For example T2K has used this 
formula in their recent publication[12] because eq.(l.l) exhibits excellent accuracy near the 
solar resonance 

So (1.3) is most likely not the true bound of validity. We would like to know to what 
extent the formula is accurate or valid. The main goal of this paper, is to mathematically 
demonstrate that there is no lower bound of A for the domain of validity. We will provide 
explicit errors of the formula, among which the main error related to the matter effect is 
only 0{si^aAA‘^). This implies that the formula is still accurate when A is close to a and 
one may apply (1.1) below the bound. 

Note that a higher order calculation in the original perturbative approach will not work 
since the series in ajA can not converge at the resonance if the branch cut singularity is 
not treated carefully. Actually a higher order correction to the formula (1.1) is computed 
in ref. [13] but the correction blows up when taking the vacuum limit A —)• 0. Thus it can 
not give a correct estimation when A is small. This is due to a lack of careful treatment of 
the branch cut singularity related to the solar resonance. 

Branch cuts in the oscillation system with the matter effect are essentially related to 
level crossings [14, 15], but less noticed before. Note that the three eigenvalues of the 
oscillation system come from the same cubic equation but they are different. The difference 
originates from the different branches in the square roots and cubic roots in the general 
solutions of a cubic equation. At a level crossing two of the eigenvalues are very close 
to each other which makes the problem quite non-perturbative and this just corresponds 
to the starting point of the branch cuts, which are called branch cut singularities. The 
branch cut singularities are essentially origins of all non-perturbativities in the oscillation 
system. In this paper, we will remove the singularity corresponding to the solar resonance 
in our analytic calculation by transformation of the eigenvalues to some singularity-free 
variables and compute the S-matrix using the Cayley-Hamilton theorem. In this way the 

^Note that for T2K, the energy range is 0.1-1.2 GeV and the spectrum peaks at 0.6 GeV[5]. A part of 
the current measured range 0.1-0.34 GeV is below the bound (1.3) which would lead to sin 20)^ > 1 in the 
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expansion. 



conventional formula will be proven to be accurate below the bound (1.3). The relation 
between the branch cut singularities and level crossings will be discussed in detail and thus 
improve our understanding of the matter effect in neutrino oscillation [16-2 9]. 

As a byproduct of our analysis, a new approximate formula is derived in this paper, 
with better accuracy. Though the exact form is a little more complicated than (1.1), for 
practical use in neutrino simulation, it is useful and covers most aspects. This is important 
considering that simulation of LBL experiments and performing require a fast and 

simple method to compute a large number of oscillation probabilities. Therefore even 
though the numerical calculation is always viable, there are still many studies on analytic 
approximation formulae for neutrino oscillation in matter [13, 30-41]. 

This paper is organized as follows. In section 2, we reformulate the original derivation 
of the formula (1.1) and numerically show the accuracy of the a-expansion in the case of 
T2K. We will see that the a-expansion for some effective parameters is actually invalid 
below 0.34GeV in T2K while the final result of the probability is very accurate. Then in 
section 3 from the viewpoint of singularities, we show that non-differentiable singularities in 
many parameters originate from the branch cut and result in the failure of the a-expansion. 
In section 4 we solve the problem rigorously and then compute the analytical error of (1.1). 
Based on the calculation in section 4, we also propose an alternative to the conventional 
formula. Their accuracies are numerically verified, which will be shown in section 5. Finally 
we conclude in section 6. 


2 The a-expansion and the accidental accuracy 


In this section, we first introduce analytic diagonalization of the 3x3 effective Hamiltonian, 
which has early been done by Zaglauer and Schwarzer [42] without any approximation. Then 
we show the a-expansion of the result from Freund’s calculation [11] and compare the ap¬ 
proximate result with the exact result (though complicated but numerically programmable) 
to see how much it deviates from the exact result. We will show that the a-expansion result 
of effective neutrino parameters are quite inaccurate and even invalid near the solar reso¬ 
nance but the final result (i.e. the assembled oscillation probability) from these parameters 
is very accurate. 

Neutrino oscillation in matter is subjected to the Schrodinger equation in the flavor 
space, 

i A|,,(i)) = hHL)), (2,1) 

where \i^{L)) denotes the flavor state of the evolving neutrino at a distance of L from the 
source and H is the Hamiltonian represented by the 3x3 matrix in the standard neutrino 
oscillation framework. 
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Here U and m, ’s are neutrino mixing matrix and masses in vacuum respectively. The second 
term in eq.(2.2) comes from the matter effect. Without the second term (i.e. Ng, = 0), the 
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solution of ( 2 . 1 ) is quite simple since the first term has already been diagonalized. So in 
vacuum, the transition amplitude of Va is just 

Sap = UaiU;^ + Ua2U*p2e'^^ + . (2.3) 


Here Sap is usually referred to as the S-matrix in neutrino oscillation. For neutrino oscil¬ 
lation in matter, we need to diagonalize (2.2) to obtain the effective mixing matrix U and 
the effective neutrino masses mfc, defined as 


H = —Udiag{m\, 




(2.4) 


Then U and rhk^ combined in the way similar to (2.3), gives the S-matrix in matter. 

The 3x3 matrix H can be analytically diagonalized by solving first the eigenvalues 
and then the corresponding eigenvectors, though the computation is complicated. 

The process can be a little simplified if we extract a dimensionless matrix M from 


and define 


Then Md is 


Md = 


_ mf Amii 

2E 2E 

(2.5) 

Md = U^MU. 

( 2 . 6 ) 

/o \ 


a + Au^.u, 

(2.7) 


\ 1 / 

where u = {Uei, Ue 2 , Ue'i) is the first row of U and is real by proper rephasing U. The cubic 
equation for the eigenvalues of Md is 


+ bX^ + cX + d = 0, 


( 2 . 8 ) 


with 

b = —1 — A — a; c = A — Au^ -I- a -I- Aa — Au\a] d = —Aau\. (2-9) 

The eigenvalues of Md (Note that M has the same eigenvalues as Md) solved from eq.(2.8) 
are 

Afc+i = -\{b + As + As), (2.10) 

where = 0 , 1 , 2 and 

1 

Ao = 6 ^ - 3c; Ai = 2b^ - 96c 27d; A 3 = _ ( 2 .II) 

Then the effective neutrino masses defined in eq.(2.4) are given by 

ml = ml + AmliAfc, 

which can be expressed in terms of a and A explicitly according to 

( 2 . 11 ). 
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( 2 . 12 ) 

eqs.(2.9), (2.10) and 






Then U can be computed by solving the corresponding eigenvectors of Afc. The reader 
may refer to [42] for the full form of eigenvectors. Once U is computed, we can extract 
effective mixing angles from the standard parametrization of U. All effective parameters 
(masses and mixing angles) expanded in a are given below[ll]: 


Ai — -{I + A — C) + 


a{C + 1 — A cos 2613)3 


12 


2C 




A 2 — CkC^2 T ^{^'^)i 

A 3 = 1(1 + .4 + C) + °(C-1+^^OOs2^33).?3 ^ 


(2.13) 

(2.14) 

(2.15) 

(2.16) 

The effective mixing angles in matter are (we use a superscript to distinguish them from 
vacuum parameters) 


where 
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In vacuum it is straightforward to get the expansion 


pvac _ 4 ,s\ 3 (^ 3 s \3 sin^ A + cos(A + h) sin A + 4si2Ci2C23a^A^ 
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(2.17) 

(2.18) 

(2.19) 

( 2 . 20 ) 

( 2 . 21 ) 


So one can replace the corresponding vacuum parameters in (2.21) with the effective pa¬ 
rameters in matter given above. This will produce the conventional formula (1.1). 

Note that the above a-expansion of effective parameters requires not only a <C 1 but 
also ajA <C 1. If we look at the effective mixing angles in (2.17-2.20), we hnd that the 
a-expansion of sin 20]^ may have a problem at ajA ~ 1. In (2.18) we see sin202^ ~ ajA 
which implies the correction from a is amplihed by 1/A, so the expansion may be not 
accurate if A is small. We compare it with the exact value from [42] in hgure 1, where 
the energy range is 0.1 — 1.2GeV and matter density is 1.3g/cm^ (the case of the T2K 
experiment). From hgure 1 we see the expansion formula hts the exact solution well only 
at E > 0.5 GeV but deviates from it quickly when E < 0.5GeV. More seriously, when the 
energy goes below 0.35 GeV then sin20]^ will be larger than 1 (the gray region). This is 
because the unitarity of U is badly violated. 

The other parameters suffering from the same problem are Ai and A2. We plot them 
with the exact solutions in hgure 2. We see that below 0.3 GeV the effective mass square 
difference Am^i = rn^{X 2 — Ai) from the exact solution (solid curves) can be several times 
that of the a-expansion (dashed curves), which also implies invalidity of the a-expansion 
at low energies. 
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Figure 1. Compare the approximate formula of sin205[^ with the exact solution. This figure shows 
the invalidity of the a-expansion of sin 20^1^ in the case of T2K. When E < 0.5GeV, it becomes 
inaccurate and for E < 0.35 GeV the result is completely invalid since the sine value should not be 
larger than l(the gray region). 



E(GeV) 


Figure 2. Compare the approximate formula of the eigenvalues Ai and A 2 with the exact solution. 

But interestingly, the formula of oscillation probability assembled from these inaccurate 
and even invalid pieces is very accurate, as shown in figure 3 where we use the same energy 
range and matter density as figure 1 and figure 2. 

One argument might be that, the accuracy of P is because sm.26'^ does not appear 
at the leading order (LO) of (1.1), but only at the next-to-leading order (NLO) and the 
next-to-next-to-leading order (NNLO) which are of order and a^, respectively. This 
suppresses the effect of the inaccuracy from sin 20^ shown in figure 1. But in figure 3 
we see that at the second and third peaks (from right to left), the NLO and NNLO are 
comparable to the LO so the accuracy can not be explained by the NLO suppression. 

It might be expected that the calculation at a higher order of a can explain this 
by finding some cancellations between errors. However, at a higher order, the accuracy 
in figures 1,2 turns out to be improved very little. Actually, as will be discussed in the 
next section, there is an underlying problem that some variables in the calculation are not 
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Figure 3. The plot shows that the conventional formula (blue solid curve) given by (1.1) is very 
accurate when used in T2K, in contrast with figure 1 and figure 2 where those effective parameters 
used to compute the probability are very inaccurate. All three order contributions (a°, and a^) 
are also plotted (dashed curves) to show that all of them are indispensable to make (1.1) accurate 
in T2K. 

differentiable at a = A = 0. For these variables, the expansion series including ajA can 
not even converge if ajA > 1. That is why a higher order calculation can not solve the 
accuracy problem. 

3 Non-differentiabilities, singularities and branch cuts in the oscillating 
system 

To reveal the key problem in the expansion, we start from an oversimplified but heuristic 
problem, series expansion of the following function 

g{a, A) = \/a'^ + A'^. (3.1) 

If a is small, but A is relatively not, then we can expand it in a, 

g(a,A) = A+^ + aO{^^). (3.2) 

Here we see ajA <C 1 is necessary to make eq.(3.2) accurate. If A is much smaller than 
a, then we can expand it in A as g{oL^ A) = a + A^ j(2a) + • • •. But what if A is close to 
a? One may think that if A is close to a, then both are small so we can make a double 
expansion of the function, 

g(a^ A) = cq T c\a T C 2 A T C)(a , aA^ A ), (3.3) 

where cq = g(0,0), ci and C 2 are the partial derivatives dag and dAg at a = A = 0. 
However, we cannot expand in this way because the partial derivatives ci and C 2 

do not exist (as one can check explicitly). Geometrically this is easy to understand since 
g = is a cone in the a — A — g space. Expansion at the tip of the cone will 

certainly fail. 
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1.5 


Figure 4. The three eigenvalues given by eq.(2.10) 
effect parameter is A = 0.1, corresponding to E = 
A^O. 


as a function of a. In the left plot, the matter 
= 1.2GeV. In the right plot we take the limit 


Actually the function (3.1) does exist in the eigenvalues of the oscillation system^ so 
in the original derivation of (1.1) ajA <C 1 lias to be assumed. If we use formulae derived 
from such an expansion but simply ignore the condition aj A <C 1, then we are at the risk 
of getting total wrong results, such as sin20(j^>l shown in Fig.l. 

So basically the question is why for the oscillation probablity this problematic expansion 
works very well. This will be answered next by branch cut singularities. 

First we look at the functions Afc(a) which are defined by the exact solution of eigen¬ 
values (2.10) and vary with a, as shown in figure 4. Note that we consider A’s as functions 
of a rather than A (or E\ since they are expanded with respect to a. 

The left plot in figure 4 shows that the eigenvalues can be very close to another at level 
crossings (corresponding to resonances), but they never really go across another. They turn 
to different directions at level crossings which implies the behavior near the resonance is 
quite non-perturbative. 

As we suppress A close to zero, the curvatures at those turning points in figure 4 become 
larger and larger. Finally the curvatures go to infinite when A —)• 0, which makes the curves 
turn suddenly at some points, then some singularities emerge. The right plot in figure 4 
shows the A —)• 0 limit. In this limit, the eigenvalues are continuous but not differentiable 
functions of a. 

In the left plot of figure 5, we plot Ai and A 2 as functions of a and A. It shows that 
there is a singularity (here we mean non-differentiable singularity) in the eigenvalues. The 
singularity is intrinsic and can not be removed by proper ordering of eigenvalues. So this 
implies that double expansion in a and A does not work. 

The intrinsic singularity in figure 5 is the kernel of the non-perturbativity problem in 
the oscillation system. It actually comes from a branch cut singularity. From eqs. (2.10) 
and (2.11) we see that Ai^ 2,3 can be analytically expressed in terms of b, c and d and then 

^The eigenvalues (and thus the corresponding oscillation parameters) contain more complicated square 
roots like \/fA + A^cl, — 2aAK where k ~ 1/3 (see, e.g., calculation in [32]), but the problem caused by 
a ~ A in the expansion is essentially the same as the simplified one in (3.1). 
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Figure 5. Singularity in the eigenvalues and its origin from branch cut singularity. In the left 
plot the eigenvalues Ai and A 2 from (2.10), as functions of a and A, have a singularity at (0,0). 
The singularity originates from the branch cut in (2.11). In the right plot we show the branch cut 
singularity of ±y^ where z = x + iy. 

further in terms of a, A and u according to (2.9). They look like smooth analytic functions 
everywhere but they are actually not. Note that there is a square root and a cubic root 
in (2.11). Functions like yT or have branch cut singularities which make them not 
analytic^ on the complex plane. 

In the right plot in figure 5 we show the two branches of where 2 = x + iy 

connecting with each other at the branch cut (the line y = 0 for x < 0). At x = 0 which is 
the end of the branch cut, there is a singularity. As one can check from (2.10) and (2.11), 
the branch cut singularity just corresponds to the singularity in the eigenvalues shown in 
the left plot. 

4 Solution 

Identifying that the singularity in the eigenvalues originates from the branch cut singularity 
makes a crucial step to solve the problem, since all branch cut singularities can be easily 
removed if the multi-branches collapse to one. For example, the branch cut singularity in 
^/z will disappear when it is squared, i.e. makes the two branches collapse and is 

singularity-free. Once the singularity is removed, ajA will not appear any more because 1 jA 
will be absorbed by some continuous and smooth functions which we will see below. After 

complex function giz) is analytic if and only if its Taylor series about 20 converges to the function 
in some neighborhood for every 20 in its domain. 
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the singularity problem is solved, we will mathematically show the conventional formula is 
accurate near the solar resonance. 

The solution of the singularity problem can be summarized by the three key steps 
below: 


1. Use the Cayley-Hamilton theorem[43-45] to express the S'-matrix only in terms of the 
eigenvalues Ai^ 2 , 3 - The eigenvectors are not needed. 

2. The singularity still exists in the eigenvalues but can be partially removed by the 
transformation 

A± = -(Ai±A2), (4.1) 

where it will be shown that A+ and A?, are singularity-free, though the singularity 
still exists in A_ since it is a branch cut singularity. 


3. It turns out that A_ only appears in cosine function and the following function 


fix) 


sinx 

X 


(4.2) 


where x is proportional to A_. Note that /(x) is smooth everywhere, even at x = 0. 

2 4 2 4 

Since /(x) = 1 — ^ + ^ -|- ... and cosx = I — ^ ^ are actually functions of 

x^ oc A^ which is singularity-free, the singularity is removed. 


In short, we will hrst make the S'-matrix only depend on (Ai,A 2 ,A 3 ) and then after the 
transformation from (Ai,A 2 ,A 3 ) to (A?_,A+,A 3 ), the singularity in the S-matrix will be 
explicitly removed. Next we will show the calculation in detail. 

The Cayley-Hamilton theorem is a theorem in linear algebra which states that if p(A) 
is the characteristic polynomial of a matrix A [for example the left-handed side of eq.(2.8)], 
then substituting the matrix A for A in this polynomial results in the zero matrix, i.e. 
piA) = 0. Take the example of eq.(2.8), this means + bM'^ + cM -|- d = 0 or 

+ cM+ d). (4.3) 


This implies = I + M + M^/2! -|- ... can be expressed by a polynomial of M with hnite 
terms since all M” with n > 3 can be converted to linear combinations of I, M, by 
eq.(4.3). So we have 

= sol + siM + S2M^, (4.4) 

where we put a —it to be used later. The coefficient sq, si and S 2 can be determined in 
various methods[43] such as the Lagrange interpolation or the Newton interpolation. They 
have been computed in [44, 45], 

So = ^[e **'''^AiA 2 (Ai — A 2 )-I- e **'’^^A 2 A 3 (A 2 — A 3 )-|-e **^^AiA 3 (A 3 — Ai)], (4.5) 

= ^[e-**"nA?-Ai) + e-**^HAi-Ai) + e-*‘^^(Ai-A?)], (4.6) 

S2 = ^[e-''^nAi-A2) + e-**^nA2-A3) + e-'‘^2(^3_^^)]^ ( 4 , 7 ) 
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where 


(4.8) 


Aa = (Ai — A2)(A2 — A3)(A3 — Ai). 


From eq.(2.1) and (2.5), the S'-matrix is 


S 


= 


. m? L 

= 


-22AM 


(4.9) 


where 



(« ) 


\ 

M = U 

1 y 

+ 

0 / 


(4.10) 


So we can identify t = 2 A to use eq.(4.4) directly. 

Now the transition amplitude of —)> I'e is Se^ = so/i 2 +siMi 2 + S 2 (M^)i 2 but because 

I is an identity matrix, can be written as two terms 


Sell = SlMu + S2(Af^)l2, 


(4.11) 


which implies we do not need to compute sq for the appearance probability. 

The denominator Aa in si and S 2 is possible to be zero, but we will show next that si 
and S 2 are not divergent at Aa = 0 and smooth (differentiable) everywhere. For example 
the singularity from Ai — A 2 = 0 can be removed by the transformation A± = ^(Ai i A 2 ). 
This singularity is the only one that confronts us in the energy range of current experiments. 
According to Vieta’s formulas for a cubic equation. 


Ai + A2 + A3 = —b, A1A2A3 = —d, (4-12) 

we have Ai + A 2 = —— A 3 and A 1 A 2 = —dAg Therefore 

A+ = li-b- A 3 ), A2_ = + dX^\ (4.13) 


where b, d are apparently free from the singularity [as shown in eq.(2.9)] and A 3 is also 
singularity-free (shown in figure 4, see also the proof in the Appendix). So A+ and A^ are 
singularity-free. Note that, however, A_ has a singularity originating from the branch cut 
singularity, which can be seen from figure 5. 

After the transformation X± = ^(Ai ± A 2 ), and S 2 are given by 


—2A+e **'^3 _|_ g **'^+[2A+ cos(A_t) -|- it{X\ -I- A?_ — A|)/(A_t)] 
A| - 2A+A3 - dAg ^ 


(4.14) 


e _|_ g **'''+[—cos(A_t)-|-it(A 3 — A+)/(A_t)] 

“ Ai - 2 A+A 3 - dAg ^ 

We see they depend only on A 3 , A+ , A_ and d which are all continuous and smooth 
functions except for A_. But since A_ only appears in cos(A_t) and /(A_t) which are 
actually functions of A^ ( note that cos(x) = 1 — ^ + ^ + ■■■ and f{x) = 1 — ^ + •■•), 

we come to the conclusion that si and S 2 are continuous and smooth functions of a and A. 


(4.15) 
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Now that we have expressed the S'-matrix in terms of (A 3 ,A+,A?_) and thus removed 
the singularity, expansion in a will not suffer from any problems. The ajA appears in 
section 2 will not appear any more if we use eqs.(4.14-4.15) to compute the probability. 
Define 

p = UesU*^, q = Ue2U*2^ (4-16) 

we have 

Sefj. = p[si + S 2 (l + A)] + qa[si + S 2 {a + A)]. (4.17) 

We call the two terms in eq.(4.17) as p term and q term respectively. From eqs.(4.14-4.15) 
we have 


p term = 


p 


Ag — 2A-|_A3 — ^^Ag 


-1 


e **'^®(A 3 — a) — e cos(A_t)(A 3 — a) 


+ite **'''+ /(A_t) (A+A3 + q;A+— d — q;A3) 


(4.18) 


q term = 


qa 


,-i 


_ 1) _ g-*tA+ cos(A_t)(A3 - 1) 


Ag — 2 A-|_A 3 — dAg 

- 2A+A3 - dAg 1 + (A3 - 1)(A3 - A+)) 


(4.19) 


So far we have not taken any approximation. Then we will use the approximation 

A 3 = 1 + 0(4^), (4.20) 

which is derived in the Appendix. With this approximation, from eq.(4.13) we have 


2A-|_ — A O' 0{si^A), A_ — A_|_ — j 


. 421 ) 

+ 0(44' ^ 

Since the p term and q term have been expressed in terms of singularity-free quantities A 3 , 
A+ and A?_, we can use (4.20,4.21) to compute them. The calculation is straightforward 
(see the Appendix) and the result is 


p term = p- 


g-2iA _ g-2iAA 


1-A 


+ 0(A4A) + 0(A2si3aA), 


term = —2iqae _|_ 0 {lS.as\^A), 


where 


A = \J {A + aY — 4AQ;Cg24’ 


(4.22) 

(4.23) 

(4.24) 


The oscillation probability is 
\p term -|- q termp 

L .ri 


= p(A) + OisfsAA) + 0(4aAA2), 


(4.25) 
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where is defined as 


p(^) = \p 


g-2fiA _ ^-2iAA 


1 - A 


A 


(4.26) 


To derive (1.1), we need to expand the modulus squared of (4.26), 

p(^) = \p termp + 2Re[p term x q term] + \q termp, (4.27) 

where \p termp, \q termp and the cross term are computed in the appendix and the result 


IS 


Ip termp = (4.28) 

\q termp = 40^5^2^12^23“”^^^ “I" O(a^A^) + OlAa^A"^) + ©(a^sisA^)), (4.29) 

2Re[p term x q term] = cos(A + sin(l — A)A ^ _|_ 0 {si5a^A). 

ss A 1 — A 


(4.30) 

Now the conventional formula (1.1) can be analytically justified near and below the 
solar resonance. Here we denote it as . We can see that the first and last terms of 
just correspond to eqs. (4.28,4.29) respectively and the middle term in corresponds to 
the cross term (4.30). Combine the analytic errors, we have 

p{B) _ p{A) ^ oisj^aA) + 0(5130^A^) + 0{a^AA^) + O(a^A^), (4.31) 


which implies that the conventional formula is accurate up to the O-terms above and the 
O-terms in (4.25). We see there is no ajA in all these O’s, so we draw the conclusion that 
the bound (1.3) which originally requires ajA <C 1 can be safely removed. The conventional 
formula is still accurate without this bound, as long as these O-terms are small. 

According to eq.(4.25), the error of is 5P^^'^ = 0{sf^AA) + 0{s‘l^aAA‘^). Taking 
T2K as an example, for E = 0.25GeV which is below the conventional domain of validity, 
we have A ~ 0.02 and A ~ 3.7. This gives sl^aAA"^ ~ 2 x 10“^ and sf 3 AA ~ 4 x 10“®. 
The error is very small and the dominant correction would be 0{s‘l^aAA‘^) if we want to 
improve the accuracy. Actually, if we only concern ourselves with the A > 1 region, then 
O(sl^aAA^) is larger than 0{sf^AA) since 513 ® — 5 ^ 3 . Therefore we expect that typically 
O(sf^aAA^) is the principal source of the error of P^^\ For the same set of parameter 
values, the four terms in eq.(4.31) have the values, Si 3 aA ~ 3 x 10“^, si 3 a^A^ ~ 2 x 10“^ 
and a^AA^ ~ a^A^ ~ 1 x 10“^. So the dominant errors of are 0{si^aA) and 

0 (si 3 a^A^). Note that they do not depend on A, which implies that the main source of 
inaccuracy of is not due to inadequately accounting for the matter effect contribution, 
but rather than due to an insufficient expansion of the small phase a A. Though the phase 
aA is small, terms quadratic in aA would not be negligible if we want to improve its 
accuracy. In conclusion, the dominant errors of P^^') and are given by 0 {s‘l^aAA‘^) 

and 0 {sl^aA) + C>(si 3 a^A^), respectively. 
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Figure 6. A comparison of the approximate oscillation formulae with the numerical solution to 
(2.1) in the T2K case. The left plot shows that all these approximate formulae are very accurate 
and the errors are negligible for practical use. The right plot shows that the residuals defined as 
|(5P| = |P—Pnumericail are Consistent with our analytic estimation of the error, which are represented 
by green and yellow shades for and respectively. 

5 Numerical verification 

As our study of the problem is originally motivated by the fact that T2K covers the solar 
resonance, we would like to numerically verify our analysis in that case first. The matter 
density in T2K is p = 2.6g/cm^[12] so we take the electron density to be Ne = 1.3A)4/cm^ 
under the assumption that for matter ZjA = 1/2 in average. 

Figure 6 shows that both and are accurate enough for practical use while the 
new formula p("^) has better accuracy than the conventional formula pi^i. We also plot the 
analytic errors according to (4.25) and (4.31) in the right panel of figure 6, using light green 
and yellow shades. The actual residuals defined as |(5P| = |P — T)iumerical| where Pnumerical 
is the numerical solution are well compatible with analytic estimation, which implies the 
errors are correctly estimated. Therefore figure 6 verifies both §pA,B) 

T2K case. 

Besides T2K, we also show the accuracies of these formulae in other accelerator neutrino 
experiments. The information of the baselines and neutrino energies are listed in table 1 
and for simplicity we take the same matter density as T2K for all the other experiments, 
since the neutrino beams in these experiments only go though the earth crust. We see 
again that in current or future accelerator neutrino experiments, the formulae are accurate 
enough for practical use and the errors are well described by our analytic estimation. 

Although computers are becoming more powerful, it is still desirable to have efficient 
methods of computation. For example, the y^-ht in a high dimensional parameter space 
(including both oscillation parameters and experiment parameters) is always extremely 
time-consuming. When nuisance parameters are being marginalized, the likelihood function 
has to be invoked an enormous amount of times to complete a sub-process of minimization 
(only for the frequentist treatment, the Bayesian approach usually needs much more com¬ 
putations). The package GLoBES [48, 49] which was designed for simulation of neutrino 
oscillation experiments has optimized the diagonalization procedure [50] combining the QL 
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Experiments 

L/km 

E/GeV 

A(E;/GeV) 

Refs 

MOMENT 

150 

~ 0.3 

0.024 

[46] 

T2K 

295 

0.6 (0.1 ^ 1.2) 

0.048 (0.008 ^ 0.096) 

[4, 5] 

MINOS 

735 

3(1 ^ 6) 

0.24 (0.08 ^ 0.48) 

[6] 

NOvA 

810 

~ 2 

0.16 

[7, 8] 

LBNE 

1300 

~ 2.5 

0.20 

[47] 


Table 1. Baseline lengths and neutrino energies of current and future accelerator neutrino ex¬ 
periments. For T2K and MINOS there are both peaks and ranges (in parentheses) of neutrino 
energies according to the references while for the other experiments we only show the general ener¬ 
gies. The electron density is = 1.3iVA/cm^ in our calculation, so we also show the values of A 
corresponding to the energies. 


decomposition with additional developed algorithms. This is not necessary for oscillations 
in constant density matter, where a simple analytical formula performs better (GLoBES 
allows users to replace the probability engine with a user-defined function). In this case, we 
recommend the use of instead. A simple test on Mathematica 8.0 with Intel Core 17 
CPU shows that 10^ evaluations of and cost^ 4.7 and 5.6 seconds, which implies 
can be computed at a speed not slower than P^^\ 

Finally there is one issue related to the solar resonance to be discussed. Strictly speak¬ 
ing, the matter effect can be safely regarded as a small perturbative effect only if 2 y/2GFNf,E 
is much less than [A <C a), i.e. only the region between the vacuum limit and the 

solar resonance can be regarded as the truly small-A region, where no physics can be 
changed greatly by the matter effect. Typically LBL accelerator neutrino experiments are 
in the region between the solar resonance and the atmospheric resonance which we can refer 
to as medium-A region. Note that originally A can not be treated perturbatively in the 
medium-A region[ll]. From the small-A region to the medium-A region, the solar mixing 
will experience a resonance. It is interesting that, according to the formulae we derived 
, the contribution from the matter effect passes through the resonance gradually without 
showing any resonances, despite the solar mixing being affected drastically in that region 
(note that the solar mixing has sizable contributions to these experiments). In other words, 
the region with a perturbative matter correction can be extended from the small-A region 
to the medium-A region for current LBL accelerator neutrino experiments. 

6 Conclusion 

The conventional formula obtained by an expansion in the mass hierarchy parameter a = 
Am 2 i/Am|]^ ~ 0.03 turns out to be very accurate near the solar resonance, as shown 
in figure 3 though the effective masses and effective mixing angles computed in the a- 

^For compiled languages which are used in practical simulation such as the C-based GLoBES package, 
the speed will be about several hundred times faster. But the ratio of the speeds of computing pM) 
varies little for different machines or different languages. 
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Figure 7. Similar plots to figure 6 but for some other accelerator neutrino experiments. For more 
details, see the T2K case in figure 6. 


expansion are inaccurate or even invalid at this region, as shown in figure 1 and figure 2. So 
it is interesting that the intermediate inaccuracies cancel each other out in the final result. 
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We have shown that the inaccuracies are because the expansion is too close to the 
branch cut singularity in the eigenvalues. This singularity is inherent in the eigenvalues so 
it cannot be removed by interchanging eigenvalues. But certain combinations of them such 
as their sum of the eigenvalues Ai + A 2 do not have the singularity, and the oscillation prob¬ 
ability only depends on these singularity-free combinations. By computing the probability 
in this way, we have analytically proven that the conventional formula is still accurate near 
the solar resonance. 

A new oscillation formula in (4.26) which might be practically useful is derived 
when we try to prove the accuracy of the conventional one. Both the conventional and 
the new formulae are very accurate in various accelerator neutrino experiments for baseline 
lengths varying from 150km (MOMENT) to 1300km (LBNE), as shown in figures 6 and 7. 
We have also estimated the analytic errors for these formulae. 

A Some details of analytic calculations 
A.l Simplify the p term and q term 

Here we show how to simplify the p term and q term step by step. All the approximations 
in the calculation should be analytically treated so we use 0() instead of ~. 

The first result we will derive is eq.(4.20). Note that the cubic equation (2.8) has the 
following identity 

b + c +d = Asi^{a — 1 ) — 1, (A.l) 

which provides a fast way to compute A 3 as follow. We assume A 3 = 1 -|- x with x <C 1 and 
replace the A in the cubic equation (2.8) with 1 -|- x. Then the leading order vanishes and 
the next-to-leading order(NLO) gives 

3x -|- x{2b -I- c) -|- O(x^) = As^ 3(1 — a), (A.2) 

which implies x = 0 ( 5 ^ 3 A) while the explicit form of x is not important here. 

Therefore from eq.(4.13) we have 

A f\/oI ^ 

2A+ = A + a + X, A2 = a2 - A 1 A 2 , A 1 A 2 = --(A.3) 

1 -|- X 

so the q term can be greatly simplified, 

qa[0{Ax) + — 2A+A3 — dX7^ + xiX^ — A+))] 

q term = ---^- 

Xj- 2A+A3 - dAg 1 
= —itqae~'''^^^ f {X-t) -|- 0{aAx) 

= sin(AA)/A + ©(aAs^gA), (A. 4 ) 
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where A = ^/{A + aY — 4Aacf2Cig or A^ = 4A?.. 


p term 


P 


(1 - A)(l - a) + 0(x} + O(aA) 

^—it\ 


e **(1 — a) — e **'''+cos(A_t)(l 


+ite * + f{X-t) {{1 — a){A — a)/2 + 0{aA) + 0{ax) + 0{X+x)) 
p [e“**(l — a) — cos(A_t)(l — a) + ite“**^+/(A_t)(l — a){A 


- a) + 0(x) 

- «)/ 2 ] 




+p[0{Ax) + 0{AaA)] 

cos(A_t) + f{X-t)^^ 

L — ./I ^ 

+p[0{Ax) + 0{AaA)]. (A.5) 


Since cos(A_t) and /(A_t) only depend on A^ = (a — A)^/4 — Q;yl(l — uf) we have 

cos(A-t) - itf{X-t) — - — = cos(—-— t) - itf {—-— t) — -h 0{aAA^). 

So hnally we get 

p term = ^ [e“** — + 0 (sf 3 AA) + 0 {sizaA/X^). (A.6) 

L — ^ 

A.2 Calculate 

After expanding the square in pA) ^ we see the square term of p equals to the leading term 
in As for the square term of q, since we have 


q = S12C13C12C23 + 0(si3), 


(A.7) 


so the differences of the corresponding term in pl-®) and the q square term is 


^2222 sin^(AA) 

4a S 12 C 12 C 23 --W square term) 


= 4a 




A2 


= 4a^ 


2a21(A2-A") + 0(si3A") 


= O(a^A^) + O(Aa^A^) + OiahuA^)), 

where we have used the following approximation for f(x) = ^ sinx 


I..9 


/(AA) - f(AA} = - A^A^) + 0(A^A^ A^A^^ 

6 


(A. 8 ) 
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The remaining term is the cross term, computed as follows 


2Re[p term x q term] 


2 Re[— 2 iapq 


g2jA _ g2iAA 


^-r(A+a)A Sin(AA) 


A 


2Re[4a(— + 0(s?3))h 


,*(A-aA+5) Sin[(l - A)A] sin(^A) 
1 -A A 


(A.9) 



+ 0(aAA3) + 0(AV)] 


(A.IO) 


+0(as?3A) + Oisua^A) + ^si 30 (a^A^,a^AA^) 

0 


(A.ll) 


where p = 8130138236 has been used. Combine the result from eqs.(A.8) and (A.IO), we 
have 


p(B) _ p(A) ^ o(8f3aA) + 0(8i3a^A^) + O(a^AA^) + O(a^A^). (A.12) 
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